## Function to calculate statistics used in the figures and 
## tables in the paper.

simulationSummary <- function(file_string){
  # file_string - location of the results of a battery's simulation
  require("plyr")
  load(file_string)
  num_qs <- ncol(sim_output$fixed_scores)
  num_obs <- length(sim_output$true_scores)
  
  errors <- array(numeric(),
                  dim = c(num_obs, num_qs, 3),
                  dimnames = list(1:num_obs,
                                  1:num_qs,
                                  c("fixed", "dynamic", "random")))
  
  errors[,,'fixed'] <- abs(sim_output$true_scores - sim_output$fixed_scores)
  errors[,,'dynamic'] <- abs(sim_output$true_scores - sim_output$adaptive_scores)
  errors[,,'random'] <- abs(sim_output$true_scores - sim_output$random_scores)
  
  variance <- array(numeric(),
                    dim = c(num_obs, num_qs, 3),
                    dimnames = list(1:num_obs,
                                    1:num_qs,
                                    c("fixed", "dynamic", "random")))

  variance[,,'fixed'] <- errors[,,'fixed']^2
  variance[,,'dynamic'] <- errors[,,'dynamic']^2
  variance[,,'random'] <- errors[,,'random']^2
  
  sum_output <- array(data = NA,
                  dim = c(num_qs, 6, 3),
                  dimnames = list(1:num_qs,
                                  c("mad", "mae", "rmse", "mpv", "meanSL", "medSL"),
                                  c("fixed", "dynamic", "random")))
  
  mad <- t(aaply(errors, 2:3, median, na.rm=TRUE))
  sum_output[,1,] <- t(mad)
  
  mae <- t(aaply(errors, 2:3, mean, na.rm=TRUE))
  sum_output[,2,] <- t(mae)
  
  calc_rmse <- function(x){
    sqrt(mean(x^2))
  }
  rmse <- t(aaply(errors, 2:3, calc_rmse))
  sum_output[,3,] <-t(rmse)
  
  mpv <- t(aaply(variance, 2:3, mean, na.rm=TRUE))
  sum_output[,4,] <- t(mpv)
  
  squaredLoss <- errors^2 + variance
  meanSL <- t(aaply(squaredLoss, 2:3, mean, na.rm=TRUE))
  sum_output[,5,] <- t(meanSL)
  
  medSL <- t(aaply(squaredLoss, 2:3, median, na.rm=TRUE))
  sum_output[,6,] <- t(medSL)
  
  return(sum_output)
}


wilcox_pval <- function(file_string){
    require("plyr")
    load(file_string)
    num_qs <- ncol(sim_output$fixed_scores)
    num_obs <- length(sim_output$true_scores)
    
    errors <- array(numeric(),
                    dim = c(num_obs, num_qs, 3),
                    dimnames = list(1:num_obs,
                                    1:num_qs,
                                    c("fixed", "dynamic", "random")))
    
    errors[,,'fixed'] <- abs(sim_output$true_scores - sim_output$fixed_scores)
    errors[,,'dynamic'] <- abs(sim_output$true_scores - sim_output$adaptive_scores)
    errors[,,'random'] <- abs(sim_output$true_scores - sim_output$random_scores)
    
    pval <- wilcox.test(errors[,,'fixed'],
                        errors[,,'dynamic'],
                        paired = FALSE)$p.value

    return(pval)
}
